Interaction broadening of Wannier functions and Mott transitions in atomic BEC 
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Superfluid to Mott-insulator transitions in atomic BEC in optical lattices are investigated for the 
case of number of atoms per site larger than one. To account for mean field repulsion between the 
atoms in each well, we construct an orthogonal set of Wannier functions. The resulting hopping am- 
plitude and on-site interaction may be substantially different from those calculated with single-atom 
Wannier functions. As illustrations of the approach we consider lattices of various dimensionality 
and different mean occupations. We find that in three-dimensional optical lattices the correction 
to the critical lattice depth is significant to be measured experimentally even for small number of 
atoms. Finally, we discuss validity of the single band model. 

PACS numbers: 03.75.Hh, 67.40.-w, 32.80.Pj, 39.25.+k 
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I. INTRODUCTION 

Numerous many-body phenomena have been recently 
demonstrated with Bose-Einstein condensates (BEC) in 
optical lattices 0, 0, 0- Number squeezing has been 
observed with 87 Rb atoms in a one-dimensional lattice 
of pancake-shaped wells pj, and superfluid to Mott- 
insulator transitions have been witnessed with such 
atoms in three-dimensional and one-dimensional optical 
lattices 2} . Such transitions were predicted by theoret- 
ical studies based on the Bose-Hubbard model jj] and 
by microscopic calculations of the model parameters for 
BEC in optical lattices [5|, |fj . 

Very important question is whether it is possible to 
observe superfluid to Mott-insulator transitions for the 
mean occupation number n larger or even much larger 
than one? Phenomenological single band Bose-Hubbard 
model indeed predicts such transitions. Previous calcu- 
lations of the model parameters J, hopping amplitude, 
and U, on-site interaction, were based on the lowest band 
Wannier functions for a single atom in an optical lattice. 
Repulsive interaction between the atoms for n > 1 may 
cause the wave function in each well to expand in all di- 
rections, not only affecting the on-site interaction U |jj 
but also strongly enhancing tunneling J between neigh- 
boring wells. This is especially significant in lower di- 
mensional lattices with transverse potential bigger than 
the lattice wells where large occupations can be achieved 
without substantial three-body collisional loss. In order 
to provide theoretical guidance for experimental observa- 
tion of Mott transitions in such systems, it is very impor- 
tant to obtain accurate critical parameters of the lattice 
potential for lattice occupations beyond unity. 

Here we show how to construct an orthogonal basis 
of Wannier functions with mean-field atomic interactions 
taken into account. We use it to obtain renormalized val- 
ues of parameters J and U, from which critical depth of 
the potential V c is calculated for various lattices of differ- 
ent dimensionality and mean occupation. For the cubic 
optical lattice with n = 2 or larger, our result is no- 



ticeably larger than that calculated without taking into 
account interaction. This increase is more pronounced 
for the anisotropic cases with stronger lattice potentials 
in one or two directions. For the case of one-dimensional 
lattice of pancake-shaped wells p| or two-dimensional lat- 
tice of tubes , our results are several times larger than 
critical values calculated from one-atom Wannier func- 
tions. This is in agreement with the experimental find- 
ings that much higher lattice potentials are needed to 
reach the transition point in such cases. 

Kohn developed variational approach to calculate elec- 
tronic Wannier functions in crystals We modify this 
procedure by minimizing on-site energy self-consistently 
taking into account interaction between atoms. 

In the last section we address validity of the single- 
band Bose-Hubbard model constructed with variational 
Wannier functions. The conditions for the model to be 
valid need to be modified from those for a single parti- 
cle case since the interaction between the particles alters 
the band structure substantially For the model to 
be valid two conditions have to be fulfilled: (i) when the 
number of particles in a well changes by one the varia- 
tional Wannier function should not change significantly 
and (ii) collective excitations of the atoms within each 
well should be less energetically favorable than atom hop- 
ping between the wells. 



II. BOSE-HUBBARD MODEL AND WANNIER 
FUNCTIONS 

For bosonic atoms located in the lattice potential V(r) 
and described by boson field operators ip(r), the Hamil- 
tonian field operator is 



H = I dr^(r) ( -^V 2 + V(r) )^(r) 
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where a s is the atoms' scattering length and m is the 
mass. To illustrate our methods we use as an example 
isotropic cubic lattice. We assume that the boson field 
operator may be expanded as ip(r) — J2ibiW(r — rj), 
where bi is the annihilation operator for an atom in the 
Wannier state of site . Substituting this expansion into 
the Hamiltonian we obtain a problem of lattice bosons. 
We consider the case when the number of atoms per cite 
rii fluctuates around average number n. This results in 
the standard Bose-Hubbard Hamiltonian 
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where the effective on-site repulsion U, the hopping am- 
plitude J and the on-site single-atom energy / are defined 
by 
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PF(r + a), (4) 
W(t), (5) 



where g = 4ira s h 2 /m and a is the lattice vector. On-site 
energy / is defined as 



/ = nl + U n{n - l)/2, 
with the bare on-site interaction Uq 



Uo = gJ dr\W(r) 



(6) 
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We assumed that the Wannier function does not change 
much for small fluctuations of the number of atoms. Off- 
site interactions are also neglected. 

In case of more than one atom per site the presence 
of other atoms does modify the Wannier function of 
an atom. Below we describe our strategy for its self- 
consistent calculation. We start with a trial wave func- 
tion localized in each well, g(r — r»). A Wannier func- 
tion corresponding to the lowest Bloch band may be con- 
structed according to Kohn's transformation: 



W(r) = ^2c H g(r - r; 

i 

dk e lkr * 



(8) 



(2tt)3 ^G(kJ' 
where the integral is over the first Brillouin zone and 

G(k) = £ / dr 5 (r) 5 (r-r i )co S (k-r < ). (9) 

i 

For an odd Wannier function, the cosine function should 
be replaced by the sine function. One can show that such 
Wannier functions are normalized and are orthogonal to 



each other for different wells. We vary the trial function 
to minimize the on-site energy / |lOj| . 

We note that another method to calculate the Wannier 
functions including interaction effects sclf-consistently 
may be used for small interactions. Starting with non- 
linear time-independent Gross-Pitaevskii equation 

-— V 2 Vr + -\i>{v)W{v) 

+ V(rty(r) = M (k)V(r), (10) 

one may calculate periodic Bloch states Uk(r) defined as 

^(r)=e lk - r Mk (r)/v^V. (11) 

by expanding them in Fourier series 

u k (x) = X A*e i2n ™/ a (12) 



and solving nonlinear system of equations. Then, a set 
of Wannier wave functions for the band in question is 
defined by 



W m (r-a) = L-Va^^r-a) 

BZ 

= L-V^^.We-*- 8 (13) 
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This procedure fails for large interactions because the 
bands develop loops and become not single- valued |||. 



III. SUPERFLUID TO MOTT-INSULATOR 
TRANSITIONS 

We consider three optical lattice systems which are 
relevant to experiments: (i) isotropic three-dimensional 
optical lattice, (ii) anisotropic three-dimensional lattices, 
and (iii) the situation when the lattice potential is present 
only in one or two directions and confinement in other 
directions is provided by relatively weak harmonic trap. 
Following standard practice, we will use the lattice period 
Tt/k, atomic mass m, and recoil energy E r = h 2 k 2 /2m as 
the basic units. 

Three pairs of counter-propagating laser beams with 
wavelength 27r/fc propagating along three perpendicular 
directions create potential 

V(r) = V x sm 2 {kx) + V y sin 2 (ky) + V z sin 2 (fcz). (14) 

Isotropic cubic lattice is created by the beam of equal 
intensity. In this case V x = V y = V z = Vq. 

Anisotropic cubic lattices can be created by choosing 
intensity of one or two beam to be much large than other. 
In this case V y = V z = V x > V x = V or V z = Vx > 
V y = V x = Vq. Below we study the case when hiuj_ 3> 
/x, where y, is the chemical potential of the atoms, thus 
the weak optical lattice is effectively one-dimensional or 
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two-dimensional and transverse motion is frozen to the 
ground state of the transverse confinement. 

Transverse motion can also be decoupled in the ex- 
perimentally relevant case when the lattice potential is 
present only in one or two directions and atoms are con- 
fined in other directions by relatively weak harmonic 
trap: Vr(r±) = ^moj\r 2 L . 

According to existing experiments, in our calculations 
through this work, we choose the 87 Rb atoms in F — 
2, m = 2 state with scattering length a s — 5.8 nm and 
the laser wavelength of 852 nm for the three- and two- 
dimensional lattices and 840 nm for the one-dimensional 
lattice. All numerical results are obtained using 21 lattice 
wells in each direction with periodic boundary condition. 
Convergence has been checked using 41 wells for some of 
the key results. 

In each case we calculate parameters of the Bose- 
Hubbard model based on the variational approach de- 
scribed in the previous section. The critical condition for 
superfluid to Mott-insulator transition has been found 
approximately as 

U/zJ = 2n + 1 + 2y/n(n + 1), (15) 

where z is the number of the nearest neighbor sites |ll| . 
By substituting the parameters into the critical condi- 
tion, we can map out the critical potential strength as a 
function of mean occupation. 

In the following, we report our findings for 
isotropic and anisotropic three-dimensional lattices, one- 
dimensional lattice of pancake wells, and two-dimensional 
lattice of tubes. 



A. Isotropic cubic lattice 

In the case of isotropic cubic lattice we choose 
variational trial function to be in the form g(r) = 
g(x)g(y)g(z), with g(u) = (1 + au 2 )e~ u ^ , where a and 
a are variational parameters. Then the Wannier function 
must also be of the product form W(r) = w(x)w(y)w(z), 
with the one-dimensional functions w(u) and g(u) related 
by the one-dimensional version of Kohn's transformation. 
All the three-dimensional integrals in Eq. (|2|)- (|15fl can 
then be reduced to one-dimensional ones, greatly simpli- 
fying the calculations. 

Our calculations proceed as following. For a given Vq 
and n, we start with certain initial parameters a and a to 
obtain a trial Wannier function through Kohn's transfor- 
mation and calculate the on-site energy /. The procedure 
is repeated by varying the parameters until the on-site 
energy / is minimized. The resulting variational Wan- 
nier function will depend on both n and Vq. If only the 
on-site single-atom energy / is minimized, one obtains 
the single-atom Wannier function Wo(r) which only de- 
pends on Vq . We find that interaction broadens Wannier 
functions, as a result Uq s is always larger than Uq, but 
we also notice that effective interaction U can be larger 
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FIG. 1: Dependence of various interaction parameters on 
number of atoms for V = 35E r . U and Uo are defined by 
J3J and Q respectively. The derivative in © is calculated 
by Chebyshev fitting to function /. Interaction parameter 
Uos calculated with single particle Wannier function is de- 
fined as Uos = g f dr|Wo(r)| 4 , where Wo(r) is a single-atom 
Wannier function. 
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FIG. 2: Hopping elements calculated with 

single particle Wannier function, Jos = 
JdrWo(r) [-ft 2 V 2 /2m + F(r)]lFo(r + a), and with the 
variational proceedure described in the text, J. Depth of the 
lattice is V = 35E r . 



than Uo (see Fig.^). So phase transition is more complex 
than we expected. 

Once the Wannier function is determined, we can cal- 
culate the Bosc- Hubbard parameters U and J. In Fig. [3 
we depict the ratio U/zJ (z = 6) as a function of the 
mean occupation n for several values of the potential 
strength Vo- The decreasing trend can be understood 
as following. The total interaction energy increases with 
n, making the Wannier function broader, hence the inter- 
action parameter U becomes smaller, J proportional to 
overlap between neighboring Wannier functions becomes 
larger, and as a result the ratio decreases. The intersec- 
tion with the line of critical condition (in Fig. [3] the line 
with positive slope obtained from Eq. H15|) ~) then yields 
the mean occupation for which these potentials are crit- 
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FIG. 3: The ratio U / zj versus mean occupation n calculated 
from the variational Wannier functions for isotropic cubic lat- 
tice. For each given parameter Vq, intersection with the solid 
line yields the mean occupation number for which the given 
Vo is critical - condition in Eg. 11511 . 



ical. For n— 1,2,3 and 4, we find the critical potentials 
to be V c = 11.95,14.32,16.25 and 18.15 respectively. A 
similar calculation can be done by using the a single- 
atom Wannier function. The critical potentials become 
11.85, 13.47, 14.61 and 15.43 for the first four mean occu- 
pations. For n = 1, the two results agree with each other 
within numerical uncertainty , and are also consistent 
with experimentally determined range for the critical po- 
tential [? ]. For n > 1, the mean field repulsion makes 
the critical potential noticeably higher. Starting from 
n = 3 the correction to the critical depth of the lattice 
has to be clearly observable experimentally and effects of 
interaction has to be taken into consideration. 



B. Anisotropic cubic lattices 

Our procedure can also be applied to the case of an 
anisotropic lattice. We model the system as a lower di- 
mensional problem with the reduced interaction parame- 
ter gd obtained by multiplying g by the integral of IV'xl 4 , 
where ip± is the single-atom ground state wave function 
in a well of the transverse potential |{|. In the har- 
monic approximation, the wave function can be found 
exactly, and the reduced interaction parameter is given 
by gi = 2£-y/V± for the quasi-one-dimensional lattice and 

■9 2 = S\f\^f^- f° r the quasi-two-dimensional lattice. In 
the calculations discussed below, we take V± = 80E r . 

To find the Wannier functions for the lower di- 
mensional lattices, we use these reduced interaction 
parameters in our procedure, replacing all the three- 
dimensional integrals in Eqs. (JSJ and JSJ by lower 
dimensional ones. The critical lattice potential V c 
calculated using such variational Wannier functions is 
depicted in Fig. 0] for the one- and two-dimensional 



FIG. 4: The critical lattice potential V c calculated from the 
variational and single-atom Wannier functions for anisotropic 
cubic lattices. The lines are guides to eyes. The dashed lines 
are for the quasi-one-dimensional and the solid lines are for 
quasi-two-dimensional cases. The triangles correspond to the 
variational and the circles to the single-particle calculations. 

models. For comparison, we also include results cal- 
culated using the one-atom Wannier function. The 
increase of critical potential due to mean-field repulsion 
on the Wannier functions is somewhat bigger in the 
lower dimensional cases. 



C. Lattices in one or two directions 

BECs in one-dimensional lattice of pancake-shaped 
wells and two-dimensional lattice of tube-shaped wells 
have been studied in experiments 1? ]. Because of the 
large transverse dimensions of such wells, many atoms 
can be held in a well without suffering too much three- 
atom collisional loss, opening the possibility of studying 
superfluid/Mott-insulator transition for relatively large 
n I n a theoretical investigation, Oosten et al [7| 

considered the interaction effect by using a transverse 
wave function in the Thomas- Fermi approximation with- 
out modifying the single-atom Wannier function in the 
lattice direction(s). Here we extend their work by con- 
sidering the interaction effect on the Wannier functions 
as well. 

For the pancake like BEC array, the transverse wave 
functions are approximated by the Thomas-Fermi wave 
function 4>tf (r_i_) of the BEC within the pancake plane, 
which is defined by 

\^T F (r ± )\ 2 = (ng,)- 1 ^ - V T (r x )), (16) 

for /i > Vt(tc±) = and vanishes otherwise. 

According to the experimental data, we take luj_ = 

19 X 27TS" 1 . 
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FIG. 6: Energy associated with hopping (process 1) has to 
be smaller than energy to excite the many-body state in well 
(process 2). Many-body excitation is schematically depicted 
as a single atom excitation. 
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FIG. 5: The critical lattice potential V c in dependence on 
mean occupation n calculated from the variational (triangle) 
and single-atom (circle) Wannier functions for: (a) the one- 
dimensional lattice with lo± = 2ir x 19 s _1 (dashed line) and 
uj± = 2-7TX 120 s _1 (solid line), and (b) two-dimensional lattice 
with lo±_ — 2tv x 24 s _1 . The lines are guides to eyes. 



We begin by writing the Wannier function in the form, 
W(y) — w(rL)(f>(r±), where (f> is the wave function for 
the transverse direction(s), and w is the Wannier func- 
tion in the lattice direction(s), both to be determined 
variationally by minimizing the on-site energy. The part 
of the on-site energy involving <j) is just the n-particle 
Gross-Pitaevskii energy in the transverse potential and 
with the interaction parameter g modified into gd by 
multiplying the integral of | w(r i ) | 4 . In the Thomas- 
Fermi approximation, this 'transverse energy' is given by 



2n- 



5ri-2 



— yjnmuj^gi/ir for the one-dimensional case and 
(Qmij^n 2 ^!) 1 / 3 for the two-dimensional case. 



Si 

fx = 

The total on-site energy is the sum of this 'transverse en- 
ergy' and n times of the single-atom energy of the lattice 
Wannier function: 



/ = fx + n 



J dr L w*{r L ) 



2m 



V 2 + V(r L ) 



Lattice Wannier function, obtained by the procedure of 
Kohn's transformation and minimization of the on-site 



energy, will be affected by the interaction because the 
'transverse energy' depends on it through the reduced 
interaction parameter gd- After w(tl) is determined vari- 
ationally, the Bose-Hubbard parameters J and U can be 
calculated immediately. In Fig. 7, we show the critical 
potential V c for the case of one-dimensional lattice with 
transverse trap frequency lu±/2tt — 19 s _1 and 120 s . 

For comparison, we also show the corresponding re- 
sults obtained using the single-atom Wannier function of 
the lattice and the Thomas-Fermi transverse wave func- 
tion. It is clear that V c is raised dramatically due to the 
broadening of the Wannier function. In the experiment 
of Ref. y, the magnetic trap potential is 19 s -1 . The 
transverse trap frequency is enhanced to 120 s^ 1 if the 
optical confining potential with Vq = 50E r is turned on, 
and the mean occupation number is n ~ 50. Evidence 
from Bragg interference pattern shows that the critical 
value of the lattice potential should be somewhat larger 
than AAE r . This observation is contradictory to the pre- 
diction based on the single-atom Wannier function, but is 
consistent with our result based on the variational Wan- 
nier function. 

In the case of two-dimensional lattice, our results 
for the critical lattice potential are shown in Fig. 
7(b) for uj ± /2tt = 24 s _1 which is used in [? ]. We 
predict V c ~ 33E r for n ~ 100, while the single-atom 
Wannier function yields V c ~ 27 E r . The largest lattice 
potential used in the experiment was 12 E r , so further 
experiment is needed to verify the theoretical predictions. 



IV. VALIDITY OF THE SINGLE-BAND MODEL 

In this section, we discuss the conditions for the single 
band Bose Hubbard model to be valid. First, we make 
general remarks and then give quantitative examples rel- 
evant for the case of the isotropic cubic lattice. 

Assumption that the boson field operator may be ex- 
panded as ip(r) — J2i biW(r — r») requires that the Wan- 
nier functions do not change substantially when the num- 
ber of atoms in a well changes by one. A good criteria for 
this condition to be fulfilled is that interaction energy cal- 
culated with the Wannier function does not change much 
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We plot the criteria from Eq. QlSj l for isotropic lattices 
on Fig. It is much smaller than unity. To estimate 
the effect of many-body excitation within a single well, 
we neglect hopping amplitude J, since close to Mott- 
insulator transition it is much smaller than atom's inter- 
action. Also for the experimentally relevant region of the 
potential depths the potential can be well approximated 
by a harmonic potential. In the harmonic potential the 
lowest many-body excited mode is associated with the 
center of mass motion - Kohn mode 0. As a result 
A ~ tkj. Since we neglect the tunneling we may start 
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FIG. 7: Relative change in the interaction energy as number 
of atoms changes by one determined by the change of the 
Wannier wave function. 



when number of particles changes by one 



\Un-Un+l\ 
U n + U n . 



< 1. 



+ 1 



(18) 



Note that the value of U can still be quite different from 
the one calculated with a single particle Wannier func- 
tion. 

When the condition is fulfilled, the second condition is 
that the excitations within the ansatz have to be the least 
energetical. That is the hopping of the atoms from well 
has to be more energetically favorable than excitation of 
atoms in each well to the many-body excited state (see 
Fig- EJ) ■ If we consider two neighboring wells the energy 
of the ground state is 

E Q = 2nI + U n(n-l). (19) 

The energy associated with hopping is 

(n-2)(n-l) + n(n + l) , s 

AEx = 2nl + £7 - t± ^ i '- -E Q = U{20) 

It has to be much smaller than the energy of the first 
excited many-body state that we denote A 



< 0.04 

3 



0.02 



U < A. 



(21) 



10 



FIG. 8: Ratio of the hopping energy to energy required to 
excite atoms in each well to the lowest many-body excited 
state. 



directly with variational form for the Wannier function 
in a well. We take W(x, y, z) = W(x)W(y)W(z), where 
W(u) = C(l + /3u 2 )e~ 7tl . Similar to previous section for 
a fixed Vq and n we minimize on-site energy /. From the 
results shown in Fig. [S] it is clear that the single-band 
model is applicable in this case: the energy associated 
with atom's hopping is much smaller than the energy 
required to excite the atoms inside of the wells. 
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